ARCH/GARCH
1.异方差问题
如果随机误差序列的方差会随着时间的变化而变化,这种情况被称为异方差。
\[Var(\epsilon_t)=h(t)\]
异方差的影响:
异方差虽然不会影响回归系数最小二乘估计的无偏性,但会影响回归系数估计的标准差和置信区间。 忽视异方差的存在会导致残差的方差会被严重低估,使得参数的显著性检验失去意义,最终导致模型的拟合精度受影响。
异方差的直观判断
残差图(方差齐性残差图,递增型残差图)
残差平方图:残差序列的方差实际上是它平方的期望 ,所以考察残差序列是 否满足方差齐性,主要考察平方序列是否平稳。
已知异方差形式:使用方差齐性变换 对于标准差与均值成正比关系的异方差序列,对数变换可以有效地实现方差齐性。 图示出现异方差现象,先尝试log看是否能够解决。
未知异方差形式:建立异方差模型。
2. ARCH模型
2.1 ARCH模型的思想
基本思想是在以前信息集下,在某一时刻一个噪声的发生是服从正态分布。该正态分布的均值为零,方差为一个随着时间变化而改变的量,即为条件异方差,并且这个随时间变化的方差是过去有限项噪声值平方的线性组合。这样就构成了自回归条件异方差模型。
2.2 ARCH模型的结构
\[ ARCH(q)\left\{ \begin{aligned} x & = \beta^T(1,x_{t-1},x_{t-2},...)+\epsilon_t \\ \epsilon_t & = \sqrt{h_t}e_t, e_t \sim N(0,1) \\ h_t & = a_0+a_1\epsilon_{t-1}^2+...+a_q\epsilon_{t-q}^2 \end{aligned} \right. \] 应当满足两个约束条件
\[ \left\{ \begin{aligned} a_0>0,\\ 0 <= a_i <1,i =1,2,...,q \\ \sum_{i=1}^qa_i<1 \end{aligned} \right. \]
- 条件一:方差必须非负,即\(Var(\epsilon_t|\epsilon_{t-1},\epsilon_{t-2},...)>=0\) 则\(a_0+a_1\epsilon_{t-1}^2+...+a_q\epsilon_{t-q}^2 >=0\)故每个参数要求非负。
- 条件二:方差平稳
2.3 ARCH模型的优缺点
优点 突破了传统时间序列模型中同方差的假设;将条件方差表达成过去干扰项的回归函数, 能够反映集聚效应;序列存在ARCH效应的时候,直接使用OLS估计会产生偏差,使用 模型能够一定程度上避免此偏差,提高参数的估计精度。
缺点 在实际应用中为了得到更好的结果,ARCH模型的阶数往往很大,参数过多,估计困 难;条件方法假设为线性函数,在现实中线性情况可能只是特例。
3.GARCH模型
3.1 模型的思想
将条件方差本身的滞后值加入\(h_t\)。即ARCH模型实际上就是\(\epsilon_t^2,\epsilon_{t-1}^2,...,\epsilon_{t-q}^2\)的q阶移动平均的MA(q)模型。而GARCH模型实际上就是\(h_t\) 关于\(h_{t-1},...,h_{t-p}\)的p阶自相关,关于\(\epsilon_t^2,\epsilon_{t-1}^2,...,\epsilon_{t-q}^2\)的q阶移动平均的ARMA(p,q)模型。
3.2 模型结构:
\[ GARCH\left\{ \begin{aligned} x_t & = f(t,x_{t-1},x_{t-2},...)+\epsilon_t \\ \epsilon_t & = \sqrt{h_t}e_t, e_t \sim N(0,1) \\ h_t & = a_0+\sum_{j=1}^p\eta_jh_{t-j}+\sum_{i=1}^q\lambda_i\epsilon_{t-i}^2 \end{aligned} \right. \] 参数满足下面的约束条件: \[ \left\{ \begin{aligned} 0<= \lambda_i<1, i=1,2,...,q\\ 0 <= \eta_j <1,j =1,2,...,q \\ 0 <= \sum_{j=1}^p\eta_j+\sum_{j=1}^q\lambda_j<1 \end{aligned} \right. \] 注: (1)\(e_t\)为原序列提取确定信息和条件异方差信息之后的残差波动,所以\(e_t\)应该是真正的 白噪声序列。\(e_t\)还有一个重要的作用,就是根据\(e_t\)的特征确定序列的分布类型。通常 假定\(e_t\)服从正态分布,如果对\(e_t\)的分布检验显示,残差序列显著拒绝正态分布假定,就需要 根据\(e_t\)的分布特征尝试其他分布。
一个完整的条件异方差模型由:均值模型,条件异方差模型和分布假定三块组成。
4.建模流程
同时关注序列信息和序列波动,首先提取序列的均值信息,然后分析残差序列中蕴含的 波动信息
- 图示法
观察适合何种模型建模,是否具有异方差性的可能。
- 线性时序建模(ARIMA, 回归拟合+残差自相关)
构建水平模型\(x_t= f(t,x_{t-1},x_{t-2},...)+\epsilon_t\) ,提取序列均值中蕴涵的相关信息. 常用ARIMA模型。建模前需要先对序列的平稳性进行检验,PP检验适用于异方差场合,故怀疑异方差的场合要重视PP检验。
建立水平模型后,检验模型残差是否通过白噪声检验,若通过,查看模型残差平方是否通过白噪声检验,若不通过,说明有异方差性,考虑继续建立ARCH/GARCH; 若一开始的残差就不通过白噪声检验,改进线性模型。
- 异方差自相关检验(LM test)
- 模型定阶
- ARCH: 看acf,pacf
- GARCH: 通常尝试(1,1),(1,2),(2,1),(2,2)
- 参数估计
- 模型检验
参数显著性检验:GARCH模型拟合完成后,首先对参数是否显著大于零进行检验,参 数显著性检验和ARMA模型相同,构建T分布统计量,在给定显著性水平下,拒绝原假 设表明该参数显著非零。
模型显著性检验:条件异方差模型拟合的效果取决于它是否将残差序列中蕴涵的异方差 信息充分提取出来。利用拟合模型估计出来的异方差 ,对残差序列和残差平方序列分 别进行标准化变换,
分布检验:\(e_t\)序列的正态性检验 在构造GARCH模型的时候,如果不特殊指定,通常默认该序列服从正态分布。故此检 验的内容之一就是检验这个分布假定是否正确。
- 模型预测
4.案例1973-2009Intel公司股票的对数收益率。
step1.时序图
0均值,无季节,趋势,所以直接拟合残差即可。
library(readr)
library(xts)
d.intel <- read_table2(
"~/Downloads/ftsdata/m-intcsp7309.txt",
col_types=cols(
.default=col_double(),
date=col_date("%Y%m%d")
))
xts.intel <- xts(
log(1 + d.intel[["intc"]]), d.intel[["date"]]
)
tclass(xts.intel) <- "yearmon"
ts.intel <- ts(c(coredata(xts.intel)), start=c(1973,1), frequency=12)
at <- ts.intel - mean(ts.intel)
ts.plot(at) 序列的均值模型是常数均值,减去均值后的残差的平方的ACF:
2.ARCH效应检验
- Ljung-Box白噪声检验
对\(a_t\),\(a_t^2\)做白噪声检验
##
## Box-Ljung test
##
## data: at
## X-squared = 18.676, df = 12, p-value = 0.09665
##
## Box-Ljung test
##
## data: at^2
## X-squared = 92.939, df = 12, p-value = 1.332e-14
统计检验方法
Portmantea Q检验
Engle拉格朗日乘子法检验(LM检验)
"archTest" <- function(x, m=10){
# Perform Lagrange Multiplier Test for ARCH effect of a time series
# x: time series, residual of mean equation
# m: selected AR order
y <- (x - mean(x))^2
T <- length(x)
atsq <- y[(m+1):T]
xmat <- matrix(0, T-m, m)
for (j in 1:m){
xmat[,j] <- y[(m+1-j):(T-j)]
}
lmres <- lm(atsq ~ xmat)
summary(lmres)
}
archTest(at,m=12)##
## Call:
## lm(formula = atsq ~ xmat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.07440 -0.01153 -0.00658 0.00395 0.35255
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.005977 0.002249 2.658 0.008162 **
## xmat1 0.093817 0.048147 1.949 0.052013 .
## xmat2 0.153085 0.048102 3.183 0.001569 **
## xmat3 0.146087 0.048614 3.005 0.002815 **
## xmat4 0.023539 0.049126 0.479 0.632075
## xmat5 0.007347 0.049107 0.150 0.881139
## xmat6 0.010342 0.047027 0.220 0.826050
## xmat7 0.057183 0.047027 1.216 0.224681
## xmat8 0.014320 0.047079 0.304 0.761149
## xmat9 0.007157 0.046968 0.152 0.878965
## xmat10 -0.019742 0.046566 -0.424 0.671810
## xmat11 -0.057537 0.046041 -1.250 0.212116
## xmat12 0.161945 0.045965 3.523 0.000473 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03365 on 419 degrees of freedom
## Multiple R-squared: 0.1248, Adjusted R-squared: 0.0997
## F-statistic: 4.978 on 12 and 419 DF, p-value: 9.742e-08
检验的p值为 9.742e-08, 高度显著,说明有ARCH效应。
3.模型建立
ARCH模型可以用acf,pacf初步判断阶数
残差平方的PACF在滞后12处较高, 另外在滞后1到3较高。 考虑建立ARCH(3)作为波动率方程。
设\(r_t\)为收益率,拟建立如下的均值方程和波动率方程: \[r_t=\mu+a_t,a_t=\epsilon_t\sigma_t, \epsilon_t i.i.d \sim N(0,1)\]
\[\sigma_t^2=\alpha_0+\alpha_1a_{t-1}^2+\alpha_2a_{t-2}^2+\alpha_3a_{t-3}^2\] 通过构造残差平方序列的自回归模型来拟合异方差函数\(\sigma_t^2\)
使用fGarch包的garchFit()函数建立ARCH模型。
library(fGarch)
mod1 <- garchFit( ~ 1 + garch(3,0), data=c(ts.intel), trace=FALSE)
#其中1表示均值方程是一个常数,
#输出结果中mu表示均值方程的均值,omega表示alpha0,alpha1为alpha1
summary(mod1)##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~1 + garch(3, 0), data = c(ts.intel), trace = FALSE)
##
## Mean and Variance Equation:
## data ~ 1 + garch(3, 0)
## <environment: 0x7ff74ae972b0>
## [data = c(ts.intel)]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 alpha2 alpha3
## 0.012567 0.010421 0.232889 0.075069 0.051994
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.012567 0.005515 2.279 0.0227 *
## omega 0.010421 0.001238 8.418 <2e-16 ***
## alpha1 0.232889 0.111541 2.088 0.0368 *
## alpha2 0.075069 0.047305 1.587 0.1125
## alpha3 0.051994 0.045139 1.152 0.2494
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 303.9607 normalized: 0.6845963
##
## Description:
## Mon Mar 8 21:14:33 2021 by user:
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 203.362 0
## Shapiro-Wilk Test R W 0.9635971 4.898647e-09
## Ljung-Box Test R Q(10) 9.260782 0.5075463
## Ljung-Box Test R Q(15) 19.36748 0.1975619
## Ljung-Box Test R Q(20) 20.46983 0.4289059
## Ljung-Box Test R^2 Q(10) 7.322136 0.6947234
## Ljung-Box Test R^2 Q(15) 27.41532 0.02552908
## Ljung-Box Test R^2 Q(20) 28.15113 0.1058698
## LM Arch Test R TR^2 25.23347 0.01375447
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -1.346670 -1.300546 -1.346920 -1.328481
得到的均值方程和波动率方程为: \[r_t=0.0126+a_t,a_t=\epsilon_t\sigma_t\] \[\sigma_t^2=0.010421+0.232889a_{t-1}^2+0.075069a_{t-2}^2+0.051994a_{t-3}^2\] 因为结果中\(\alpha_2\)和\(\alpha_3\)的估计值是不显著的,可以拟合ARCH(1)模型为波动率方程
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~1 + garch(1, 0), data = c(ts.intel), trace = FALSE)
##
## Mean and Variance Equation:
## data ~ 1 + garch(1, 0)
## <environment: 0x7ff74b9e22c0>
## [data = c(ts.intel)]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1
## 0.013130 0.011046 0.374976
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.013130 0.005318 2.469 0.01355 *
## omega 0.011046 0.001196 9.238 < 2e-16 ***
## alpha1 0.374976 0.112620 3.330 0.00087 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 299.9247 normalized: 0.675506
##
## Description:
## Mon Mar 8 21:14:33 2021 by user:
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 144.3783 0
## Shapiro-Wilk Test R W 0.9678175 2.670321e-08
## Ljung-Box Test R Q(10) 12.12248 0.2769429
## Ljung-Box Test R Q(15) 22.30705 0.1000019
## Ljung-Box Test R Q(20) 24.33412 0.2281016
## Ljung-Box Test R^2 Q(10) 16.57807 0.08423723
## Ljung-Box Test R^2 Q(15) 37.44349 0.001089733
## Ljung-Box Test R^2 Q(20) 38.81395 0.007031558
## LM Arch Test R TR^2 27.32897 0.006926821
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -1.337499 -1.309824 -1.337589 -1.326585
这个模型的aic比上个模型差一些。
为了进行模型验证,可以计算标准化残差\(\tilde \alpha_t = \frac{a_t}{\sigma_t}\) 计算标准化残差\(\tilde \alpha_t\)
resi <- residuals(mod2, standardize=TRUE)
#标准化残差的时序图
plot(ts(resi, start=start(ts.intel), frequency=frequency(ts.intel)),
xlab="年", ylab="标准化残差")检验标准化残差的分布。
##
## Shapiro-Wilk normality test
##
## data: resi
## W = 0.96782, p-value = 2.67e-08
另:fGarch包对建模见过自带若干诊断图
获得拟合的均值用fitter()函数,
## 1 2 3 4 5 6
## 0.009999835 -0.150012753 0.067064079 0.082948635 -0.110348491 0.125162849
用predict函数做超前若干步的预测
## meanForecast meanError standardDeviation
## 1 0.01313003 0.1090513 0.1090513
## 2 0.01313003 0.1245215 0.1245215
## 3 0.01313003 0.1298482 0.1298482
## 4 0.01313003 0.1317901 0.1317901
## 5 0.01313003 0.1325109 0.1325109
## 6 0.01313003 0.1327802 0.1327802
## 7 0.01313003 0.1328811 0.1328811
## 8 0.01313003 0.1329188 0.1329188
## 9 0.01313003 0.1329330 0.1329330
## 10 0.01313003 0.1329383 0.1329383
## 11 0.01313003 0.1329403 0.1329403
## 12 0.01313003 0.1329411 0.1329411
预测包括均值的预测(显然是用\(\mu\)预测的)、 均值预测的标准误差、 波动率的预测(是用的值滚动计算的)。 波动率长期预测接近于ARCH模型的无条件标准差。
GARCH模型
ARCH模型用来描述波动率能得到很好的效果, 但实际建模时可能需要较高的阶数, 比如§17.5.3的欧元汇率波动率建模用了11阶的ARCH模型。
GARCH模型的定阶方法研究不多, 一般用试错法尝试较低阶的GARCH模型, 如GARCH(1,1), GARCH(2,1), GARCH(1,2)等。 许多情况下GARCH(1,1)就能解决问题。
继续用上面的数据建模
library(fGarch, quietly = TRUE)
mod1 <- garchFit(~ 1 + garch(1,1), data=ts.intel, trace=FALSE)
summary(mod1)##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~1 + garch(1, 1), data = ts.intel, trace = FALSE)
##
## Mean and Variance Equation:
## data ~ 1 + garch(1, 1)
## <environment: 0x7ff74a45c030>
## [data = ts.intel]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 beta1
## 0.01126568 0.00091902 0.08643831 0.85258554
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.0112657 0.0053931 2.089 0.03672 *
## omega 0.0009190 0.0003888 2.364 0.01808 *
## alpha1 0.0864383 0.0265439 3.256 0.00113 **
## beta1 0.8525855 0.0394322 21.622 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 312.3307 normalized: 0.7034475
##
## Description:
## Mon Mar 8 21:14:34 2021 by user:
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 174.904 0
## Shapiro-Wilk Test R W 0.9709615 1.030282e-07
## Ljung-Box Test R Q(10) 8.016844 0.6271916
## Ljung-Box Test R Q(15) 15.5006 0.4159946
## Ljung-Box Test R Q(20) 16.41549 0.6905368
## Ljung-Box Test R^2 Q(10) 0.8746345 0.9999072
## Ljung-Box Test R^2 Q(15) 11.35935 0.7267295
## Ljung-Box Test R^2 Q(20) 12.55994 0.8954573
## LM Arch Test R TR^2 10.51401 0.5709617
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -1.388877 -1.351978 -1.389037 -1.374326
对标准化残差及其平方的白噪声检验都通过了。而且AIC值更好。条件分布的正太性检验仍然不通过。 模型可以写成: \[r_t = 0.0113 + a_t, a_t = \sigma_t\epsilon_t, \epsilon_t i.i.d. \sim N(0,1)\]
\[\sigma_t^2 = 0.00092+0.0864a_{t-1}^2+0.852\sigma_{t-1}^2\]
如果均值非常数,需要回归拟合,使用ARMA+GARCH模型
这里没有对arma后的白噪声做过多的要求,直接根据aic判断
library(quantmod)
library(lattice)
library(timeSeries)
library(rugarch)
loadSymbols("^GSPC", from="1950-01-01")
spReturns = diff(log(Cl(GSPC)))
spReturns[as.character(head(index(Cl(GSPC)),1))] = 0
# Create the forecasts vector to store the predictions
windowLength = 500
foreLength = length(spReturns) - windowLength
forecasts <- vector(mode="character", length=foreLength)
for (d in 0:foreLength) {
# Obtain the S&P500 rolling window for this day
spReturnsOffset = spReturns[(1+d):(windowLength+d)]
# Fit the ARIMA model
final.aic <- Inf
final.order <- c(0,0,0)
for (p in 0:5) for (q in 0:5) {
if ( p == 0 && q == 0) {
next
}
#通过aic寻找最优的arma阶数
arimaFit = tryCatch( arima(spReturnsOffset, order=c(p, 0, q)),
error=function( err ) FALSE,
warning=function( err ) FALSE )
if( !is.logical( arimaFit ) ) {
current.aic <- AIC(arimaFit)
if (current.aic < final.aic) {
final.aic <- current.aic
final.order <- c(p, 0, q)
final.arima <- arima(spReturnsOffset, order=final.order)
}
} else {
next
}
}
#对残差拟合garch模型
spec = ugarchspec(
variance.model=list(garchOrder=c(1,1)),
mean.model=list(armaOrder=c(final.order[1], final.order[3]), include.mean=T),
distribution.model="sged"
)
fit = tryCatch(
ugarchfit(
spec, spReturnsOffset, solver = 'hybrid'
), error=function(e) e, warning=function(w) w
)
# If the GARCH model does not converge, set the direction to "long" else
# choose the correct forecast direction based on the returns prediction
# Output the results to the screen and the forecasts vector
if(is(fit, "warning")) {
forecasts[d+1] = paste(index(spReturnsOffset[windowLength]), 1, sep=",")
print(paste(index(spReturnsOffset[windowLength]), 1, sep=","))
} else {
fore = ugarchforecast(fit, n.ahead=1)
ind = fore@forecast$seriesFor
forecasts[d+1] = paste(colnames(ind), ifelse(ind[1] < 0, -1, 1), sep=",")
print(paste(colnames(ind), ifelse(ind[1] < 0, -1, 1), sep=","))
}
}
# Output the CSV file to "forecasts.csv"
write.csv(forecasts, file="forecasts.csv", row.names=FALSE)
# Input the Python-refined CSV file
spArimaGarch = as.xts(
read.zoo(
file="forecasts_new.csv", format="%Y-%m-%d", header=F, sep=","
)
)
# Create the ARIMA+GARCH returns
spIntersect = merge( spArimaGarch[,1], spReturns, all=F )
spArimaGarchReturns = spIntersect[,1] * spIntersect[,2]
# Create the backtests for ARIMA+GARCH and Buy & Hold
spArimaGarchCurve = log( cumprod( 1 + spArimaGarchReturns ) )
spBuyHoldCurve = log( cumprod( 1 + spIntersect[,2] ) )
spCombinedCurve = merge( spArimaGarchCurve, spBuyHoldCurve, all=F )
# Plot the equity curves
xyplot(
spCombinedCurve,
superpose=T,
col=c("darkred", "darkblue"),
lwd=2,
key=list(
text=list(
c("ARIMA+GARCH", "Buy & Hold")
),
lines=list(
lwd=2, col=c("darkred", "darkblue")
)
)
)其他:https://rpubs.com/englianhu/binary-Q1FiGJRGARCH
第四章条件异方差作业
- 查找格力电器(000651)日收盘价,并尝试用ARCH模型对其波动率建模。
da <- read.csv('~/Desktop/000651.csv',fileEncoding = "GBK")
library(lubridate)
ts.stock <- zoo(da[,c('收盘价')],
as.POSIXct(da$日期))
plot(as.xts(ts.stock), type="l",
multi.panel=TRUE, theme="white",
major.ticks="years",
grid.ticks.on = "years")## Series: ts.stock
## ARIMA(1,1,2) with drift
##
## Coefficients:
## ar1 ma1 ma2 drift
## 0.7970 -1.1364 0.1616 0.0501
## s.e. 0.0455 0.0630 0.0569 0.0279
##
## sigma^2 estimated as 22.15: log likelihood=-1479.45
## AIC=2968.9 AICc=2969.02 BIC=2989.96
##
## Call:
## arima(x = ts.stock, order = c(5, 1, 3))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3
## 1.0996 -0.8540 0.4276 -0.0236 -0.0129 -1.4495 1.1917 -0.6777
## s.e. 0.3374 0.3378 0.0970 0.0808 0.0966 0.3350 0.4464 0.1544
##
## sigma^2 estimated as 21.39: log likelihood = -1472.69, aic = 2963.39
##
## Box-Pierce test
##
## data: model$residuals
## X-squared = 0.012938, df = 1, p-value = 0.9094
##
## Box-Pierce test
##
## data: model$residuals^2
## X-squared = 47.769, df = 1, p-value = 4.796e-12
ARCH模型可以用acf,pacf初步判断阶数。
使用fGarch包的garchFit()函数建立ARCH模型。
library(fGarch)
mod1 <- garchFit( ~ 1+arma(5,1,3)+ garch(2,0), data=ts.stock, trace=FALSE)
#其中1表示均值方程是一个常数,
#输出结果中mu表示均值方程的均值,omega表示alpha0,alpha1为alpha1
summary(mod1)##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~1 + arma(5, 1, 3) + garch(2, 0), data = ts.stock,
## trace = FALSE)
##
## Mean and Variance Equation:
## data ~ 1 + arma(5, 1, 3) + garch(2, 0)
## <environment: 0x7ff74c456388>
## [data = ts.stock]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu ar1 ar2 ar3 ar4 ar5 omega alpha1
## -0.62940 1.00000 -0.76323 0.70764 0.72471 -0.66544 1.82423 1.00000
## alpha2
## 1.00000
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu -0.62940 0.65074 -0.967 0.333436
## ar1 1.00000 0.07707 12.976 < 2e-16 ***
## ar2 -0.76323 0.08582 -8.894 < 2e-16 ***
## ar3 0.70764 0.09704 7.292 3.05e-13 ***
## ar4 0.72471 0.10167 7.128 1.02e-12 ***
## ar5 -0.66544 0.07975 -8.344 < 2e-16 ***
## omega 1.82423 0.52616 3.467 0.000526 ***
## alpha1 1.00000 0.08155 12.262 < 2e-16 ***
## alpha2 1.00000 0.09491 10.536 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## -1302.881 normalized: -2.605762
##
## Description:
## Mon Mar 8 21:14:35 2021 by user:
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 328093.6 0
## Shapiro-Wilk Test R W 0.5015763 0
## Ljung-Box Test R Q(10) 11.96909 0.2871305
## Ljung-Box Test R Q(15) 14.70098 0.473163
## Ljung-Box Test R Q(20) 19.84824 0.4674582
## Ljung-Box Test R^2 Q(10) 0.1092745 1
## Ljung-Box Test R^2 Q(15) 0.1653461 1
## Ljung-Box Test R^2 Q(20) 0.2268573 1
## LM Arch Test R TR^2 0.124615 1
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## 5.247523 5.323386 5.246890 5.277292
resi <- residuals(mod1, standardize=TRUE)
#标准化残差的时序图
plot(ts(resi, start=start(ts.stock), frequency=frequency(ts.stock)),
xlab="年", ylab="标准化残差")4.查找沪深300指数日数据,并尝试用ARCH模型对其波动率建模。
da <- read.csv('~/Desktop/399300.csv',fileEncoding = "GBK")
library(lubridate)
ts.stock <- zoo(da[,c('收盘价')],
as.POSIXct(da$日期))
plot(as.xts(ts.stock), type="l",
multi.panel=TRUE, theme="white",
major.ticks="years",
grid.ticks.on = "years")## Series: ts.stock
## ARIMA(3,1,1) with drift
##
## Coefficients:
## ar1 ar2 ar3 ma1 drift
## -0.7462 0.0368 0.0986 0.8037 4.7986
## s.e. 0.1151 0.0561 0.0490 0.1088 2.6901
##
## sigma^2 estimated as 2912: log likelihood=-2695.75
## AIC=5403.49 AICc=5403.66 BIC=5428.77
##
## Call:
## arima(x = ts.stock, order = c(3, 1, 1))
##
## Coefficients:
## ar1 ar2 ar3 ma1
## -0.7348 0.0481 0.1061 0.7989
## s.e. 0.1150 0.0559 0.0486 0.1089
##
## sigma^2 estimated as 2901: log likelihood = -2697.31, aic = 5404.61
##
## Box-Pierce test
##
## data: model$residuals
## X-squared = 0.013309, df = 1, p-value = 0.9082
##
## Box-Pierce test
##
## data: model$residuals^2
## X-squared = 5.5184, df = 1, p-value = 0.01882
ARCH模型可以用acf,pacf初步判断阶数。
使用fGarch包的garchFit()函数建立ARCH模型。
library(fGarch)
mod1 <- garchFit( ~ 1+arma(3,1,1)+ garch(1,0), data=ts.stock, trace=FALSE)
#其中1表示均值方程是一个常数,
#输出结果中mu表示均值方程的均值,omega表示alpha0,alpha1为alpha1
#根据系数的显著性做调整
summary(mod1)##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~1 + arma(3, 1, 1) + garch(1, 0), data = ts.stock,
## trace = FALSE)
##
## Mean and Variance Equation:
## data ~ 1 + arma(3, 1, 1) + garch(1, 0)
## <environment: 0x7ff74c0a60c0>
## [data = ts.stock]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu ar1 ar2 ar3 omega alpha1
## 13.990177 1.000000 0.010294 -0.012689 2269.843557 0.268760
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 13.99018 19.09086 0.733 0.46367
## ar1 1.00000 0.05943 16.828 < 2e-16 ***
## ar2 0.01029 0.06513 0.158 0.87441
## ar3 -0.01269 0.04427 -0.287 0.77438
## omega 2269.84356 185.78754 12.217 < 2e-16 ***
## alpha1 0.26876 0.08187 3.283 0.00103 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## -2695.309 normalized: -5.390618
##
## Description:
## Mon Mar 8 21:14:38 2021 by user:
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 207.0629 0
## Shapiro-Wilk Test R W 0.9589577 1.391679e-10
## Ljung-Box Test R Q(10) 14.35463 0.1574218
## Ljung-Box Test R Q(15) 19.44239 0.1943819
## Ljung-Box Test R Q(20) 23.17653 0.28021
## Ljung-Box Test R^2 Q(10) 29.10891 0.001196442
## Ljung-Box Test R^2 Q(15) 37.02431 0.001255208
## Ljung-Box Test R^2 Q(20) 38.56193 0.00755475
## LM Arch Test R TR^2 25.38049 0.01311926
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## 10.80524 10.85581 10.80495 10.82508
resi <- residuals(mod1, standardize=TRUE)
#标准化残差的时序图
plot(ts(resi, start=start(ts.stock), frequency=frequency(ts.stock)),
xlab="年", ylab="标准化残差")使用rugarch包进行garch,便于一起预测arma和garch
library(rugarch)
model.garch = ugarchspec(mean.model=list(armaOrder=c(3,1,1),include.mean=TRUE),
variance.model=list(model='sGARCH',garchOrder=c(1,0)))
#model Valid models (currently implemented) are “sGARCH”, “fGARCH”, “eGARCH”, “gjrGARCH”, “apARCH” and “iGARCH” and “csGARCH”.
model.garch.fit = ugarchfit(data=ts.stock, spec=model.garch )
model.garch.fit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,0)
## Mean Model : ARFIMA(3,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 3099.700409 1.983654 1562.6212 0.000000
## ar1 0.099566 0.052889 1.8826 0.059761
## ar2 0.954577 0.007328 130.2577 0.000000
## ar3 -0.047632 0.046344 -1.0278 0.304050
## ma1 1.000000 0.000025 39994.3468 0.000000
## omega 2265.939234 187.502340 12.0849 0.000000
## alpha1 0.236764 0.074350 3.1844 0.001450
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 3099.700409 2.081529 1489.1460 0.000000
## ar1 0.099566 0.042965 2.3174 0.020484
## ar2 0.954577 0.005886 162.1868 0.000000
## ar3 -0.047632 0.037819 -1.2595 0.207856
## ma1 1.000000 0.000027 37074.4144 0.000000
## omega 2265.939234 337.919059 6.7056 0.000000
## alpha1 0.236764 0.117688 2.0118 0.044241
##
## LogLikelihood : -2689.533
##
## Information Criteria
## ------------------------------------
##
## Akaike 10.786
## Bayes 10.845
## Shibata 10.786
## Hannan-Quinn 10.809
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.0134 0.9078
## Lag[2*(p+q)+(p+q)-1][11] 5.0447 0.9498
## Lag[4*(p+q)+(p+q)-1][19] 8.8548 0.6626
## d.o.f=4
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.09406 0.7591
## Lag[2*(p+q)+(p+q)-1][2] 0.42373 0.7305
## Lag[4*(p+q)+(p+q)-1][5] 4.28760 0.2203
## d.o.f=1
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[2] 0.6541 0.500 2.000 0.41865
## ARCH Lag[4] 5.0922 1.397 1.611 0.08359
## ARCH Lag[6] 9.1950 2.222 1.500 0.02158
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.7898
## Individual Statistics:
## mu 0.12571
## ar1 0.05968
## ar2 0.09902
## ar3 0.05362
## ma1 0.29883
## omega 0.30924
## alpha1 0.11060
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.69 1.9 2.35
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.2952 0.7680
## Negative Sign Bias 0.4298 0.6675
## Positive Sign Bias 0.8137 0.4162
## Joint Effect 1.0815 0.7815
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 38.80 0.004688
## 2 30 45.16 0.028388
## 3 40 60.16 0.016358
## 4 50 74.40 0.011110
##
##
## Elapsed time : 0.592417
##
## Box-Pierce test
##
## data: resi
## X-squared = 0.52614, df = 1, p-value = 0.4682
##
## Box-Pierce test
##
## data: resi^2
## X-squared = 6.954, df = 1, p-value = 0.008363
- 预测
样本外预测
##
## *------------------------------------*
## * GARCH Model Forecast *
## *------------------------------------*
## Model: sGARCH
## Horizon: 3
## Roll Steps: 0
## Out of Sample: 0
##
## 0-roll forecast [T0=2021-01-07]:
## Series Sigma
## T+1 5536 65.12
## T+2 5536 57.18
## T+3 5553 55.14
样本内预测
out_of_sample <- round(length(ts.stock)/2)
dates_out_of_sample <- tail(index(ts.stock), out_of_sample)
garch_fit <- ugarchfit(spec = model.garch, data = ts.stock, out.sample = out_of_sample)
coef(garch_fit)## mu ar1 ar2 ar3 ma1 omega
## 3.110350e+03 1.889647e-01 7.442907e-01 6.876422e-02 8.579572e-01 1.984805e+03
## alpha1
## 4.397274e-02
garch_fore <- ugarchforecast(garch_fit, n.ahead = 1, n.roll = out_of_sample-1)
forecast_value <- xts(garch_fore@forecast$seriesFor[1, ], dates_out_of_sample)
forecast_volatility <- xts(garch_fore@forecast$sigmaFor[1, ], dates_out_of_sample)
# plot of ts.stock value
plot(cbind("fitted" = fitted(garch_fit),
"forecast" = forecast_value,
"original" = ts.stock),
col = c("blue", "red", "black"), lwd = c(0.5, 0.5, 2),
main = "Forecast of ts.stock", legend.loc = "topleft") #拟合不错# plot of volatility
plot(cbind("fitted volatility" = sigma(garch_fit),
"forecast volatility" = forecast_volatility,
"log-returns" = log(ts.stock)),
col = c("blue", "red", "black"), lwd = c(2, 2, 1),
main = "Forecast of volatility of synthetic log-returns", legend.loc = "topleft")R Session Information
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.5
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] zh_CN.UTF-8/zh_CN.UTF-8/zh_CN.UTF-8/C/zh_CN.UTF-8/zh_CN.UTF-8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] rugarch_1.4-4 lubridate_1.7.9 fGarch_3042.83.2
## [4] fBasics_3042.89.1 timeSeries_3062.100 timeDate_3043.102
## [7] xts_0.12.1 zoo_1.8-8 readr_1.4.0
## [10] forecast_8.13 dplyr_1.0.2
##
## loaded via a namespace (and not attached):
## [1] mclust_5.4.7 Rcpp_1.0.6
## [3] mvtnorm_1.1-1 lattice_0.20-41
## [5] digest_0.6.27 lmtest_0.9-38
## [7] truncnorm_1.0-8 R6_2.5.0
## [9] evaluate_0.14 ggplot2_3.3.2
## [11] pillar_1.4.6 rlang_0.4.10
## [13] spatial_7.3-12 curl_4.3
## [15] fracdiff_1.5-1 TTR_0.24.2
## [17] nloptr_1.2.2.2 SkewHyperbolic_0.4-0
## [19] Matrix_1.2-18 rmarkdown_2.5
## [21] labeling_0.4.2 stringr_1.4.0
## [23] munsell_0.5.0 compiler_4.0.3
## [25] numDeriv_2016.8-1.1 xfun_0.18
## [27] DistributionUtils_0.6-0 pkgconfig_2.0.3
## [29] urca_1.3-0 htmltools_0.5.1.1
## [31] nnet_7.3-14 Rsolnp_1.16
## [33] tidyselect_1.1.0 tibble_3.0.4
## [35] bookdown_0.21 quadprog_1.5-8
## [37] crayon_1.4.1 MASS_7.3-53
## [39] GeneralizedHyperbolic_0.8-4 grid_4.0.3
## [41] nlme_3.1-149 gtable_0.3.0
## [43] lifecycle_1.0.0 magrittr_2.0.1
## [45] scales_1.1.1 rmdformats_0.3.7
## [47] KernSmooth_2.23-17 quantmod_0.4.17
## [49] stringi_1.5.3 farver_2.0.3
## [51] tseries_0.10-47 ellipsis_0.3.1
## [53] generics_0.0.2 vctrs_0.3.4
## [55] spd_2.0-1 tools_4.0.3
## [57] glue_1.4.2 purrr_0.3.4
## [59] hms_0.5.3 ks_1.11.7
## [61] yaml_2.2.1 colorspace_1.4-1
## [63] knitr_1.30